The impact of short term synaptic depression and 
stochastic vesicle dynamics on neuronal variability 
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Abstract Neuronal variability plays a central role in 
neural coding and impacts the dynamics of neuronal 
networks. Unreliability of synaptic transmission is a 
major source of neural variability: synaptic neurotrans- 
mitter vesicles are released probabilistically in response 
to presynaptic action potentials and are recovered stochas- 
tically in time. The dynamics of this process of vesicle 
release and recovery interacts with variability in the 
arrival times of presynaptic spikes to shape the vari- 
ability of the postsynaptic response. We use continuous 
time Markov chain methods to analyze a model of short 
term synaptic depression with stochastic vesicle dynam- 
ics coupled with three different models of presynaptic 
spiking: one model in which the timing of presynaptic 
action potentials are modeled as a Poisson process, one 
in which action potentials occur more regularly than a 
Poisson process and one in which action potentials oc- 
cur more irregularly. We use this analysis to investigate 
how variability in a presynaptic spike train is trans- 
formed by short term depression and stochastic vesicle 
dynamics to determine the variability of the postsynap- 
tic response. We find that regular presynaptic spiking 
increases the average rate at which vesicles are released, 
that the number of vesicles released over a time win- 
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dow is more variable for smaller time windows than 
larger time windows and that fast presynaptic spiking 
gives rise to Poisson-like variability of the postsynaptic 
response even when presynaptic spike times are non- 
Poisson. Our results complement and extend previously 
reported theoretical results and provide possible expla- 
nations for some trends observed in recorded data. 
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1 Introduction 



Variability of neural activity plays an important role in 
population coding and network dynamics 21 . Random 
fluctuations in the number of action potentials emitted 
by a population of neurons affects the firing rate of 
downstream cells [501H5] . In addition, spike count vari- 
ability over both short and long timescales can impact 
the reliability of a rate-coded signal [TC] . It is therefore 
important to understand how this variability is shaped 
by synaptic and neuronal dynamics. 

Several studies examine the question of how intrin- 
sic neuronal dynamics interact with variability in presy- 
naptic spike timing to determine the statistics of a post- 
synaptic neuron's spiking response, but many of these 
studies do not account for dynamics and variability in- 
troduced at the synaptic level by short term synaptic 
depression and stochastic vesicle dynamics. Synapses 
release neurotransmitter vesicles probabilistically in re- 
sponse to presynaptic spikes and recover released vesi- 
cles stochastically over a timescale of several hundred 
milliseconds [551152]. The dynamics and variability in- 
troduced by short term depression and stochastic vesi- 
cle dynamics alter the response properties of a postsy- 
naptic neuron ^^^^^^^b^EJl^M^!^ 
and therefore play an important role in information 
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transfer [5MIMM1IMH5] , neural coding [jnEMHZZMIl 
and network dynamics [^|^|71|f^|^|H] . Under- 
standing how variability in presynaptic spike times in- 
teract with short term depression and stochastic vesicle 
dynamics to determine the statistics of the postsynaptic 
response is therefore an important goal. 

In this study, we use a model of short term synap- 
tic depression with stochastic vesicle dynamics to ex- 
amine how variability in a presynaptic input is trans- 
ferred to variability in the synaptic response it pro- 
duces. We use the theory of continuous-time Markov 
chains to construct exact analytical methods for calcu- 
lating the the statistics of the postsynaptic response to 
three different presynaptic spiking models: one model 
with Poisson spike arrival times, one with more regular 
spike arrival times, and one with more irregular spike 
arrival times. We find that depressing synapses shape 
the timescale over which neuronal variability occurs: 
the number of neurotransmitter vesicles released over 
a time interval is highly variable for shorter time win- 
dows, but less variable for longer time windows when 
variability is quantified using Fano factors. Addition- 
ally, we find that when presynaptic inputs are highly ir- 
regular (Fano factor greater than 1), synaptic dynamics 
cause a reduction in Fano factor, consistent with pre- 
vious studies 26,25,18,19 . On the other hand, when 
presynaptic input is more regular (Fano factor less than 
1), synaptic dynamics often cause an increase in Fano 
factor. This observation suggests a mechanism through 
which irregular and Poisson-like variability can be sus- 
tained in spontaneously spiking neuronal networks [541 
[5Tll8ll9l [35lfTT| , which complements previously proposed 
mechanisms [58 6152,2911331. 

2 Methods 

We begin by introducing the synapse model used through- 
out this study. We then proceed by analyzing the statis- 
tics of the synaptic response to three different input 
models. 

2.1 Synapse model 

A widely used model of depressing synapses [5"71l2"ll5"5l 
I54"ll4"8] does not capture stochasticity in vesicle recovery 
and release. As a result, this model underestimates the 
variability of the synaptic response [TDTI3] . For this rea- 
son, we use a more detailed synapse model that takes 
stochastic recovery times and probabilistic release into 

account [501[n21[211H3] • 

We consider a presynaptic neuron with spike train 
I(t) — J2j S(t — tj) that makes M functional contacts 



onto a postsynaptic cell. Here, tj is the time of the 
jth presynaptic action potential. Define m(t) to be the 
number of contacts with a readily releasable neuro- 
transmitter vesicle at time t (so that < m(t) < M). 
For simplicity, we assume that each contact can release 
at most one neurotransmitter vesicle in response to a 
presynaptic spike. When a presynaptic spike arrives, 
each contact with a releasable vesicle releases its vesi- 
cle independently with probability p r . After releasing 
a vesicle, a synaptic contact enters a refractory period 
during which it is unavailable to release a vesicle again 
until it recovers by replacing the released vesicle. The 
recovery time at a single contact is modeled as a Pois- 
son process with rate l/r n . Equivalently, the duration of 
the refractory period is exponentially distributed with 
mean t u . 

Define Wj to be the number of contacts that release 
a vesicle in response to the presynaptic spike at time 
tj (so that Q < Wj < m(tj) < M where m(tj) = 
lim. t - m(t)). The synaptic response is quantified by 
the marked point process 

x{t) = ^2w j d{t-t j ). 

3 

Since the signal observed by the postsynaptic cell is de- 
termined by x(t), we quantify synaptic response statis- 
tics in terms of the statistics of x(t) in our analysis. The 
process x(t) can be convolved with a post-synaptic re- 
sponse kernel to obtain the conductance induced on the 
postsynaptic cell [23J . The effects of this convolution on 
response statistics is well understood |53j , so we do not 
consider it here. 

This model can be described more precisely using 
the equation [33] 

dm(t) = -dN x (t) + dN u (t) 

where dN u (t) — u(t)dt is the increment of an inhomo- 
geneous Poisson process with instantaneous rate that 
depends on m(t) through (dN u (t)) | m(t)) = dt(M - 
m(t))/T u (here, (• | •) denotes conditional expectation), 
N x (t) = f Q x(s)ds is the number of vesicles released up 
to time t, and each Wj is a binomial random variable 
with mean p r m(tj) and variance m(tj)p r (l — p r ). 

2.2 Statistical measures of the presynaptic spike train 
and the synaptic response 

We focus on steady state statistics in this article, and 
therefore assume that the presynaptic spike trains are 
stationary and that the synapses have reached statis- 
tical equilibrium. The intensity of a presynaptic spike 
train is quantified by the mean presynaptic firing rate, 

r in = (I(t)) = (Nj(T)/T) 
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where (•) denotes the expected value and 



2.3 Model analysis with Poisson presynaptic inputs 



N I (T) = / I(s)ds 
Jo 

represents the number of spikes in the time interval 
[0,T]. Temporal correlations in the presynaptic spike 
times are quantified by the auto-covariance, 

R iD (T)=COv{I(t),I{t + T)), 

and the variability in the presynaptic spike train is 
quantified by its Fano factor, 



var(7V/(T)) 



For much of this work, we will focus on Fano factors 
over large time windows which, through a slight abuse 
of notation, we denote by F ul = limx->oo F U1 (T). To 
compute Fano factors, we will often exploit their rela- 
tionship to auto-covariance functions, 



Fi n (T) = — / fli„(r)(l - \r\/T)dr 

r in J—T 



and 



Rin(T)dT. 



(1) 



(2) 



The statistics of the synaptic response, x(t), are de- 
fined analogously to the statistics of I(t). The steady 
state rate of vesicle release is defined as 

r x = (x(t)) = (N X {T)/T) 

where N X (T) — x(s)ds represents the number of 
vesicles released in the time interval [0,T]. Temporal 
correlations in the synaptic response are quantified by 
the auto-covariance, 

R x {t) = cov(x(t),x(t + T)) 

and response variability is quantified by the Fano factor 
of the number of vesicles released, 



F X (T) 



var(A^(T)) 



As above, we define F x — limT->oo F X (T) and note that 
F X {T) = — f R x {T){l)dr and F x = f R x (r)dT. 



(3) 



We first consider a homogeneous Poisson input, 
with rate ri n . The input auto-covariance for this model 
is given by i?i n (r) = ri n 8(r) and the Fano factor is given 
by Fi n (T) = 1 for any T > 0. The mean rate of vesicle 
release for this model is given by 

Mp r r- m 



P r rinT u + 1 



which saturates to M/t u for large presynaptic rates, r; n . 
A closed form approximations to the auto-covariance 
function of the response for this Poisson input model 
are derived in [431136] (see also [17]) and consist of a 
sum of a delta function and an exponential, 



= Dr x S(r) - Er x e-W T °, 
where the mass of the delta function is given by 
2pr (n n T u + M - 1) + 2 - p r 2 n n T u 



(4) 



D = 



> 0, (5) 



(2 - p r )p r r m T u + 2 
the timescale of the exponential decay is given by 



to 



1 + PrnnTu ' 

and the peak of the exponential is given by 
p r r in ((M - 2) Pr + 2)t„ + 2(M - l)p r + 2 



E = 



M (2 - p r )Prn n T u + 2M 



(6) 



It can easily be checked that E > whenever M > 1, 
< p r < 1, rin > 0, and t u > so that the peak of the 
exponential in ^ is negative. For finite T, the Fano 
factor, F X (T), is given by 



F X (T) = D - 2Et - [e ^ -1 
and, in the limit of large T, 



2Erl 
T 



F X =D- 2Ern 



(7) 



(8) 



To test the accuracy of these approximations, ex- 
act solutions can be found numerically using standard 
methods for the analysis of continuous-time Markov 
chains, as described for alternate input models below. 
This analysis is a special case of the analysis for the 
"regular" input model described below that is achieved 
by taking 9 = 1. Alternatively, exact numerical results 
can be achieved by taking r s = rj for the "irregular" 
input model. In figures showing results for the Poisson 
input model, we plot the closed form approximations 
described above along with exact numerical results ob- 
tained using the regular input model with 0=1. 
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2.4 Model analysis with irregular presynaptic inputs 



then calculated explicitly to obtain 



Spike trains measured in vivo often exhibit irregular 
spiking statistics indicated by Fano factors larger than 
1 [4lH5l[3l[TTj . To describe the synaptic response to ir- 
regular inputs, we use a model of presynaptic spiking in 
which the instantaneous rate of the presynaptic spike 
train, I(t), randomly switches between two values, r s 
and 77 > r s , representing a slow spiking state and a 
fast spiking state. The time spent in the slow state 
before transitioning to the fast state is exponentially 
distributed with mean t s . Likewise, the amount of time 
spent in the fast state before switching to the slow state 
is exponentially distributed with mean 77. Transition 
times are independent from one another and from the 
spiking activity. Between transitions, spikes occur as a 
Poisson process. 

To find ri n , i?i n (r), and F- m , we represent this model 
as a doubly stochastic Poisson process. Define r(t) 6 
{r s ,rf} to be the instantaneous firing rate at time t. 
Then r(t) is a continuous time Markov chain |31j on 
the state space r — (?* s , 77) with infinitesimal generator 
matrix 

~-1/t, 1/ts ' 

Clearly, r(t) spends a proportion t s /(t s + 77) of its 
time in the slow state (defined by r(t) = r s ) and a pro- 
portion Tf/(r s + Tf) of its time in the fast state (defined 
by r(t) = 77). This gives a steady-state mean firing rate 
of 

_ r s T s + 7777 

r in - 

Ts+Tf 

At non-zero lags (r ^ 0), the auto-covariance of a 
doubly stochastic Poisson process is the same as the 
auto-covariance of r(t) [43J, which we can compute us- 
ing techniques for analyzing continuous time Markov 
chains. For r > 0, we have 

{r(t)r(t + r)) = r s Pr(r(i) = r s )(r(t + r) | r{t) = r s ) 
+ 77 Pr(r(t) = r f )(r(t + r)\r(t)=r f ) 



Pr(r(i + t) = r s \r(t) = r s ) 



A 



r s r s 

T s +Tf 
r f T f 



(r(t + T)\r(t) =r s ) (9) 
(r(t + r)\r(t)=r f ) 



T s + Tf 

where (■)■) denotes conditional expectation and 

(r(t + t) I r(t) = r s ) = r s Pr(r(t + t) = r s | r(t) = r s ) 
+ r f (l- Pr(r(t + t) = r s \ r(t) = r s )). 

The probability in this expression can be written in 
terms of an exponential of the generator matrix, A, and 



77 e T f Ts +t s 



Tf+T s 

where [v]k denotes the kth component of a vector, v. 
An identical calculation can be performed to obtain an 
analogous expression for (r(t + t) \ r(t) = r f) . Combin- 
ing these with Eq. © gives 



(r(t)r(t + T)) 



_ TfT s (r f -r s ) -^-^ 

(Tf+T s f 



rl . 



For positive t, we have Ri n (T) — (r(t)r(t + T)) —rf n . As 
with all stationary point processes Ri n {T) — R ln (— t) 
and i?i n (r) has a Dirac delta function with mass ri n 
at the origin [Tl]. Thus, the auto-covariance of I(t) is 
given by 



Rin(r) = r in S(T) + 



TfT„(rf -r 3 ) 2 ^-M_M 

{Tf+Tsf 



(10) 



For finite T, the Fano factor, F; n (T), can be computed 
using Eqs. (JXJ) and (fTUj) . In the limit of large T, we can 
use Eqs. (0) and (fTUj) to obtain a closed form expression, 



Fin = 1 



2r / V A 2 (r / -r s ) 2 

{Tf+Tsf {TfTf +r s T s ) 



(11) 



Poisson spiking is recovered by setting rf = r Sl 
Tf = 0, or r s = 0. For any other parameter values (i.e., 
when rf ^ r s and 77, t s > 0), it follows from Eq. 
that Fi n (T) > 1 for any T. Therefore this input model, 
hereafter referred to as the "irregular spiking" model, 
represents spiking that is more irregular than a Poisson 
process. 

The analysis in [33] used to derive closed form ex- 
pressions for the response statistics with Poisson in- 
puts cannot easily be generalized to derive expressions 
with non-Poisson inputs like those considered here. In- 
stead, we analyze the synaptic response for the irregu- 
lar input model using techniques for analyzing contin- 
uous time Markov chains. First note that the process 
b(t) — {m(t), r(t)) is a continuous-time Markov chain on 
the discrete state space {0, 1, ... , M} x {r s , 77}. Here, 
m{t) denotes the size of readily releasable pool and r(t) 
represents the instantaneous presynaptic rate (which 
switches between r a and rf). We enumerate all 2(M+1) 
elements of this state space and denote the jth ele- 



ment of this enumeration as Fj 
1,...,2(M+1). 



( m j": r i) for J 
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The infinitesimal generator, B, of b(t) is a 2(M + 
1) x 2(M + 1) matrix with off-diagonal terms defined 
by the instantaneous transition rates, 

B j>k = lim 1 Pr(6(t + ft) = A | 6(t) = A), (12) 

and with diagonal terms chosen so that the rows sum 
to zero: B 



3 -J 



To fill the matrix B, we consider each type of tran- 
sition that the process b(t) undergoes. Vesicle recovery 
events occur at the instantaneous rate (M — m(t))/T u 
and increment the value of unit) by one vesicle. There- 
fore 

1 M — m 

lim -Pr(b{t + h) = (m+l,r)\b(t) = (m,r)) = 

h->-0 ft t u 

for to G {0, . . . , M — 1} and r G {r s , r/}. Vesicle release 
events occur at the instantaneous rate r(t) and decre- 
ment the value of m(t) by a random amount k with a 
binomial distribution so that 

lim i Pr(&(t + ft) = (to - fc, r) I &(t) = (to, r)) = 
/t-^oft 



(TO-fc)! 



i—k 



for to £ {1, ... , M}, k G {0, . . . , to}, and r G {r s , r/}. 
The value of r(t) switches from r s to r/ with instanta- 
neous rate l/r s so that 

lim - Pr(6(i + ft) = (m, r/) | 6(t) = (to, r s )) = — 

fl— >-0 ft ' T s 

and, similarly, 

lim ~ Pr(6(t + ft) = (to, r s ) | b(t) = (m, r/)) = — . 

/i— >0 ft T / 

These four transition types account for all of the tran- 
sitions that bit) undergoes. They can be used to fill 
the off-diagonal terms of the matrix B. The diagonal 
terms are then filled to make the rows sum to zero, as 
discussed above. 

Once a the infinitesimal generator matrix, B, is ob- 
tained, the probability distribution of b(t) given an ini- 
tial distribution p(0) is given by 



Pit) 



tB 



p(0). 



The stationary distribution, poi of b(t) is given by the 
vector in the one-dimensional null space of B with ele- 
ments that sum to one [oTj . 

The instantaneous rate of vesicle release, conditioned 
on the current state of r(t) and m(t), is given by 

(x(t) | r(t) — r, m(t) = m) = rp r m. 



Averaging over r and to in the steady state gives 

) 

\pa]jrjp r mj 



2(M+1) 

E 

3=1 



where [•] j denotes the jth element. The auto-covariance, 
R x (r), has a Dirac delta function at r = 0. We separate 
this delta function from the continuous part by writing 
Rx{t) = A x 8{t) + R£(t) where R£(t) is a continuous 
function. The area of the delta function can be found by 
conditioning on the current state of r(t) in the steady 
state to get 

A x = lim (x(t) 2 dt) 



lim (dN 2 x (t)/dt) 



2{M+l) 

E 

3 = 1 



r iNj l im (wl\m(t k ) = rnj) 



(13) 



where w k is the number of vesicles released by the 
kth presynaptic spike. Conditioned on the size, m(t^), 
of the readily releasable pool immediately before the 
presynaptic spike arrives, Wk has a binomial distribu- 
tion with second moment, 

lim (w k | m(t7) — raj) — m,jp r (l — p r ) + m 2 p 2 

which can be substituted into Eq. to calculate A x . 

All that remains is to calculate the continuous part, 
R£(t), of R x (t). First note that, for r > 0, 

lim (x(t)x(t + r)) 

t— >oo 

2(M+1) 

= \po] i Pv(b(t + T) = r j \b(t) = r i ) (14) 

x (dN x (t)dN x (t + r) | b(t) = r h b{t + r)= r 3 )/dt 2 . 

The second term in Eq. ([T4l can be computed as 

Pr(6(< + r) = r, | b{t) = A) = K^eJ, = [e^fc, , 

where is the 2(M+ 1) x 1 vector whose zth element is 
1 and all other elements are zero, which represents an 
initial distribution concentrated at Pj. The last term in 
Eq. (|T4"]) is given by 

(dN x it)dN x it + t) i bit) = r iy b(t +r) = r^/dt 2 = 



nrjPr mjTOj. 

Finally, R x ir) = lim t _ 
so that 



,(x(t)x(t + \t\)) - r 2 x for r / 



2(Af+l) 



TiTjPr milTlj 



(i 
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which can be computed efficiently using matrix multi- 
plication. The response Fano factor, F x , can then be 
found by integrating R x {t) according to Eqs. ^ and 



2.5 Model analysis with regular presynaptic inputs 

We now consider a spiking model that gives Fano fac- 
tors smaller than 1 and therefore spike trains that are 
more regular than Poisson processes. We achieve this by 
defining a renewal process with gamma-distributed in- 
terspike intervals (ISIs). Such a process can be obtained 
by first generating a Poisson process, J^k ^ ~ s fe) with 
rate r = 9r- ln for some positive integer 9, then keep- 
ing only every 9th spike to build the spike train I(t). 
More precisely, the first spike of the gamma process is 
obtained by choosing an integer, k, uniformly from the 
set and defining defining t\ = Sk- The re- 

maining spikes are defined by tj + \ — sjg + k to obtain 
the stationary renewal process, I(t) — J2j $(t — tj) P2]- 
Clearly, this process has rate r; n since the original 
Poisson process has rate #n n and a proportion 1/9 of 
these spikes appear in I(t). The auto-covariance is given 

by m 



Rin(r) = r in S(T) + r h 



Jk(r) 



(15) 



where 



fk(t) 



t ke-u er .)ke e -9r in t 



(k9 - 1)1 

is the density of the waiting time between the first spike 
and the (fc + l)st spike (i.e., the duration of k consecu- 
tive ISIs). 

For finite T, the Fano factor, F; n (T), can be com- 
puted using Eqs. ([T]) and (fT5j) . In the limit of large T, 
we can use Eq. ([2]) or use the fact that, for renewal pro- 
cesses, F in = var(ISI)/(ISI} 2 where var(ISI) = l/(rf n 0) 
is the variance and (ISI) = l/r- m is the mean of the 
gamma distributed ISIs [T3]. This gives 

p. I 
111 ~ 9- 

Poisson spiking is recovered by setting 9 = 1. When 
9 > 1, we have that F in (T) < 1 for any T. Therefore 
this model, hereafter referred to as the "regular" input 
model, represents spiking that is more regular than a 
Poisson process. 

The synaptic response with the regular input model 
can be analyzed using methods similar to those used for 
the irregular model. We introduce an auxiliary process, 
q(t), that transitions sequentially through the state space 
{!,..., 6}. Once reaching 9, q(t) transitions back to 



state 1. Transitions occur as a Poisson process with 
rate 9r ln . The waiting times between transitions from 
q = 8 to q = 1 are gamma distributed. Thus, to recover 
the regular input model, we specify that each transi- 
tion from 1 = 8 to q = 1 represents a single presynaptic 
spike. The process g(t) = (m(t),q(t)) is then a con- 
tinuous time Markov chain on the discrete state space 
{1,...,6>} x {0,...,M}. We enumerate all 9(M + 1) 
elements of this space and denote the jth element as 
r j = (m j ,q j ) tor j = 1,...,9(M+1). 

The infinitesimal generator, G, which is a 9(M + 
1) x 9(M + 1) matrix is defined analogously to the ma- 
trix B in Eq. (|T2l) above. The elements of G can be 
filled using the following transition probabilities. As for 
the irregular input model, vesicle recovery occurs as a 
Poisson process with rate (M — m(t))/T u so that 

1 M — 777 
lim - Pi(g(t+h) = (to+1, q) \ g(t) = (m, <?)) = 

for 777 = 0, . . . , M and q = 1, . . . , 9. Transitions that 
increment q{t) occur with instantaneous rate, 9r- m so 
that 

lim i Pi{g{t + h) = (m,q+l)\ g(t) = (m, qj) = 6r m 

for q = 1, . . . , 9 — 1 and m = 0, . . . , M . The only other 
transitions are those from q(t) = 9 to q(t) = 1, which 
represent a presynaptic spike and are therefore accom- 
panied by a release of vesicles. The transitions con- 
tribute the following, 



lim - Pr(<?(t + h) 
h^o h 



9n 



(l,m-k)\g(t) 

Pr k (l- Pr ) m - k 



(9,m)) 



{m-k)V 

for 777 G {0, 1, ... , M} and k £ {0, . . . , m}. These tran- 
sition rates can be used to fill the off-diagonal terms 
of the matrix G. The diagonal terms are then filled so 
that the rows sum to zero. The stationary distribution, 
Po, of g(t) = (m(t),q(t)) is given by the vector in the 
one-dimensional null space of G with elements that sum 
to one. 

A proportion [po^fc) of time is spent in state (m(t), q(t)) 
(k,6) where 7(fc) represents the index of the element 
(k,9) in the enumeration chosen for r (i.e., the index, 
j, at which Fj = (k,9)). In that state, the transition 
to q(t) = 1 occurs with instantaneous rate 9r m and re- 
leases average of p r m(t) vesicles. Thus, the mean rate 
of vesicle release is given by 



r x = 2_ J r inPrk [p ] 



fe=i 



As above, we separate the auto-covariance into a 
delta function and a continuous part by writing R x (r) = 
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A x 8{t) + R x (t) where R^{t) is a continuous function. 
The area of the delta function at the origin is given by 



M 



A X =2_ J # r in[Po] 7 0) ( fe Fr(l - Pr) + k 2 p r 2 ) 



k=l 



by an argument identical to that used for the irregular 
input model above. Also by a similar argument used for 
the irregular input model, we have that 



A x S{t) - t\ + e 2 rf n X>o] 7 (fc) 



M 



k,l=l 



7(0,7(fc)' 



Pr kl. 



2.6 Parameters used in figures 

Theoretical results are obtained for arbitrary parameter 
values, but for all figures we use a set of parameter val- 
ues that are consistent with experimental studies. For 
synaptic parameters, we use t u — 700 ms and p r = 0.5 
consistent with measurements of short term depression 
in pyramidal-to-pyramidal synapses in the rat neocor- 
tex [571122] . We also choose M = 5 which is within the 
range observed in several cortical areas [6]. 

The Poisson presynaptic input model is determined 
completely by its firing rate and the regular input model 
is determined completely by its firing rate and Fano 
factor. Presynaptic firing rates and Fano factors are 
reported on the axes or captions of each figure. The 
irregular input model has four parameters that deter- 
mine the firing rate and Fano factor. In all figures, we 
set Tb = t s = 1.315/c, r\, — 37c, and r s — 3c which 
gives a Fano factor of F[ n = 20.0017 « 20 for any value 
of c (from Eq. (JTTJ)) . Changing c effectively scales the 
timescale of presynaptic spiking, hence scaling r- m , with- 
out changing F m . 



3 Results 

We analyze the synaptic response to different patterns 
of presynaptic inputs using a stochastic model of short 
term synaptic depression in which a presynaptic neuron 
makes M functional contacts onto a postsynaptic neu- 
ron [50,22,24JT5]. The input to the presynaptic neuron 
is a spike train denoted by I(t). Neurotransmitter vesi- 
cles are released probabilistically in response to each 
presynaptic spike. Specifically, a contact with a readily 
available vesicle releases this vesicle with probability p r 
in response to a single presynaptic spike. After a synap- 
tic contact has released its neurotransmitter vesicle, it 
enters a refractory state where it is unable to release 
again until the vesicle is replaced. The duration of this 
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Fig. 1 Rate of vesicle release as a function of presynaptic 
firing rate for various presynaptic Fano factors. The rate 
of vesicle release, r x , is an increasing function of presynap- 
tic firing rate, ri n . Vesicle release is slower for the irregular 
spiking model than for the Poisson and regular spiking model. 



refractory period is an exponentially distributed ran- 
dom variable with mean t u , so that vesicle recovery is 
Poisson in nature. 

We are interested in how the statistics of the presy- 
naptic spike train determine the statistics of the synap- 
tic response. The presynaptic statistics are quantified 
using the presynaptic firing rate, rj n , the presynaptic 
auto-covariance function, i?j n (r), and the Fano factor, 
Fi n (T), of the number of presynaptic spikes during a 
window of length T. Similarly, we quantify the statistics 
of the synaptic response using the mean rate of vesicle 
release, r x , the auto-covariance of vesicle release, R x (t), 
and the Fano factor, F X (T), of the number of vesicles 
released during a window of length T. We will espe- 
cially focus on Fano factors over large time windows 
and define F in = lim T ^oo F in (T), F x = lim T ^oo F X (T) 
accordingly. See Methods for more details. 

We begin by considering the effect of F{ B on the 
mean rate of vesicle release, r x . We then examine the 
dependence of F X (T) on the length, T, of the time win- 
dow over which vesicle release events are counted. Fi- 
nally, we show that short term synaptic depression pro- 
motes Poisson-like responses to non-Poisson presynap- 
tic inputs. 



3.1 Irregular presynaptic spiking reduces the rate at 
which neurotransmitter vesicles are released 

We first briefly investigate the dependence of the rate 
of vesicle release, r x , on the rate and variability of the 
presynaptic spike train, as measured by r- ln and Fj n re- 
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Fig. 2 Synaptic response to an irregular and a regular presynaptic spike train. A) An irregular spike train, I(t), drives 
a depressing synapse. Each vesicle release event is indicated by a vertical bar with height indicating the number of vesicles 
released (here, all events release just one vesicle). Each time a vesicle is released, the number of available vesicles, m(t), is 
decremented accordingly. Vesicle recovery increments m(t) and occurs randomly in time (vesicle recovery events indicated by 
filled triangles). B) Same as (A) except for a more regular presynaptic spike train. Note that, even though the same number of 
spikes occur in (A) and (B), the irregular spike train is less effective in releasing vesicles. This occurs because all vesicles are 
depleted by the first few spikes in a burst and subsequent spikes in that burst are unable to release vesicles. For illustrative 
purposes, we set M = 3 in this figure. 



spectively. Vesicle release rate generally increases with 
n n , but saturates to r x = M/t u whenever p r r ln 3> 1/t u 
since synapses are depleted in this regime (Fig. [T}. 

When presynaptic spike times occur as a Poisson 
process (so that F; n = 1), the mean rate of vesicle 
release is given by r x = Mp r ri n /(p r n n Tu + 1) [2"2"UT§1 
Interestingly, vesicle release is slower for more ir- 
regular presynaptic spiking and faster for more regular 
presynaptic spiking even when presynaptic spikes ar- 
rive at the same mean rate (Fig. [TJ also see [IB])- This 
can be understood by noting that, for the irregular in- 
put model, spikes arrive in bursts of higher firing rate 
followed by durations of lower firing rate. Vesicles are 
depleted by the first few spikes in a burst and subse- 
quent spikes in that burst arc ineffective and therefore 
essentially "wasted" spikes (Fig. [3]A) . When presynap- 
tic spikes arrive more regularly, more vesicles are re- 
leased on average (Fig. [2j3) . 



3.2 Variability in the number of vesicles released in a 
time window decreases with window size 

We now consider how the the variability of the synap- 
tic response to a presynaptic input depends on the 
timescale over which this variability is measured. We 
quantify the variability of the synaptic response us- 
ing the Fano factor, F X (T), which is defined to be the 
variance-to-mean ratio of the number of vesicles re- 
leased in a time window of length T (see Methods) 



and can be calculated from an integral of the auto- 
covariance function, R x {t), using Eq. (|3]). 

The auto-covariance of a Poisson presynaptic spike 
train is simply a delta function at the origin, i?j n (r) = 
fmS{T), and the Fano factor over any window size is 
therefore equal to one, F in (T) = 1 (Figs. [3p and HP). 
The auto-covariance of the synaptic response when presy- 
naptic inputs are Poisson consists of a delta function 
at the origin surrounded by a double-sided exponential 
with a negative peak (see Eq. Q and Fig. |3p) that de- 
cays with timescale tq — r u j (I + p r ri n T u ). The fact that 
the auto-covariance is negative away from r = implies 
that the Fano factor, F X (T), is monotonically decreas- 
ing in the window size, T (see Eq. Q and Fig. HP). 
For small T, the mass of the delta function at the ori- 
gin dominates the integral in Eq. ([3]) so that the Fano 
factor is approximately equal to the ratio of this mass 
to the mean rate, r x , at which vesicles are released. As 
T increases, the negative mass of the exponential peak 
subtracts from the positive contribution of the delta 
function and decreases the Fano factor. In particular, 
F X {T) w D - FT + 0(T 2 ) where Dr x is the mass of 
the delta function in R x (t) and —Er x is the peak of 
the exponential in R x {t) (see Eqs. ([5]) and ([5])). As T 
continues to increase, F X (T) monotonically decreases 
towards its limit, F x := liniT->oo Fx(T) = D — IEtq. 
Thus, short term synaptic depression converts a Fano 
factor that is constant with respect to window size into 
one that decreases with window size (Fig. HP,D). 
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Fig. 3 Presynaptic auto-covariance functions and auto-covariance of vesicle release for three input models. Auto- 
covariance functions of presynaptic spike trains (top row) and the synaptic response they evoke (bottom row) for three 
different presynaptic input models: A,B)irregular (_F; n = 20), C,D) Poisson (F; n = 1) and E,F) regular (F- ln = 0.1). Each 
auto-covariance function has a Dirac delta function at the origin that is not depicted here. Dotted line in (D) is from the 
approximation in Eq. Q and solid line is from exact calculation obtained using numerics for the regular input model with 
9=1 (see Methods), but the two are nearly indistinguishable. i?i n (r) has units (spikes/sec) 2 and R x (r) has units (vesicles/sec) 2 . 
Short term depression introduces negative temporal correlations even when presynaptic spike trains are temporally uncorrelated 
(C,D) or positively correlated (A,B). 



When presynaptic spike times are not Poisson, the 
statistics of the postsynaptic response cannot be de- 
rived analytically using the methods utilized for the 
Poisson input model. Instead, we use the fact that the 
synapse model can be represented using a continuous 
time Markov chain, which can be analyzed to derive 
expressions for the response statistics in terms of an 
infinitesimal generator matrix (see Methods). 

Irregular presynaptic spiking (i.e., inputs with F m > 
1) is achieved by varying the rate of presynaptic spik- 
ing randomly in time to produce a doubly stochas- 
tic Poisson process (see Methods). For this model, the 
input auto-covariance is a delta function at the ori- 
gin surrounded by an exponential peak (see Eq. ITUl 
and Fig. EJ^). The input Fano factor therefore increases 
with window size (see Eq. [TT1 and Fig. H^-). The posi- 
tive temporal correlations exhibited in the input auto- 
covariance function are canceled by the temporal de- 
correlating effects of short term synaptic depression [2lfl 
HZ5TI24) . For the parameters chosen in this study, this 
de-correlation outweighs the positive presynaptic cor- 
relations so that the auto-covariance function of the 
response is negative away from r = (Fig. [33), al- 
though parameters can also be chosen so that temporal 
correlations in the response are small and positive [25 . 
As with the Poisson input model, negative temporal 
correlations cause the response Fano factor to decrease 



with window size (Fig. HJ3). Thus short term synaptic 
depression and stochastic vesicle dynamics can convert 
a presynaptic Fano factor that increases with window 
size into one that decreases. 

Regular presynaptic spiking is achieved by generat- 
ing a renewal process with gamma-distributed intcrspike- 
intervals. The input auto-covariance function for this 
model exhibits temporal oscillations (Eq. ([T5|) and Fig.[3j3) 
and the Fano factor generally decreases with window 
size (Fig. [4(3). Perhaps unsurprisingly, the auto-covariance 
function of the synaptic response exhibits oscillations 
and the response Fano factor decreases with window 
size (Figs.[3F andiF). 

For all three input models, the variability of the 
synaptic response is larger over shorter time windows 
and smaller over larger time windows. A postsynap- 
tic neuron that is in an excitable regime will generally 
respond most effectively to inputs that exhibit more 
variability over short time windows [16,47,39,38]. In 
addition, rate coding is often more efficient when spike 
counts over larger time windows are less variable |64) . 
Thus, the dependence of F X (T) on window size is espe- 
cially efficient for the neural transmission of rate-coded 
information |24j . 

In addition to the temporal dependence of F X (T) 
introduced by short term depression, note that the re- 
sponse Fano factor for the irregular input model is sub- 
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Fig. 4 Presynaptic and response Fano factor as a function of window size for three input models. Presynaptic and 
response Fano factors, Fi n (T) and F X (T), as a function of the window size over which inputs or vesicles are counted (see 
Methods), obtained by applying Eq. ([TJ to the auto-covariance functions in Fig. [3] Short term depression causes response Fano 
factor to decrease with window size even when presynaptic Fano factor increases with window size (A,B) or is independent of 
window size (C,D). Also, response Fano factors are near 1 even when presynaptic Fano factors are not (B and F). 



stantially smaller than the input Fano factor (Fig. 
Conversely, the response Fano factor for the regular in- 
put model is larger than the input Fano factor (Fig. [4]). 
For both models, the response Fano factor is substan- 
tially nearer to 1 than the input Fano factor. We explain 
this phenomenon next. 

3.3 Depleted synapses exhibit Poisson-like variability 
even when presynaptic inputs are highly non-Poisson 

We now investigate the dependence of the variability 
in synaptic response on the rate and variability of the 
presynaptic input. Since we have already discussed the 
dependence of F X (T) on T above, we will focus here 
on the Fano factor calculated over long time windows, 
F x =lim T ^oo F X {T). 

We first consider parameter regimes where the ef- 
fective rate of presynaptic inputs is much slower than 
the rate of vesicle recovery (p r r- m <C 1/t„). In such 
a regime, each contact is likely to recover between two 
consecutive presynaptic spikes and therefore all M con- 
tacts are likely to have a vesicle ready to release when 
each spike arrives (Fig. In this limit, the number 
of vesicles released by each spike is an independent bi- 
nomial variable with mean (wj) = p r M and variance 
var(wj) = Mp r (l — p r ). The number, N X (T), of vesi- 
cles released in a time window of length T can then be 
represented as a sum of N- m (T) independent binomial 
random variables (i.e., a random sum). The mean of 



this sum is given by (N X (T)} = (N in (T))(wj}, which 
implies that r x = Mp r r m in this limit. Similarly, the 
variance of this sum is given by [31] va.r(N x (T)) — 
{N in (T)}va,r(Wj) + (Wj) 2 v&r(N in (T)), which implies 

lim F X (T) = ( Wj )F in (T) + vax( Wj ) / ( Wj ) (16) 

ri n ->0 

= l+p r (MF in (T)-l). 

Eq. (1161) is verified for the Poisson input model by tak- 
ing ri n — > in Eqs. (J7J. For the irregular and regular 
input models, Eq. (|16[) should be interpreted heuristi- 
cally, as it was derived heuristically. A counterexample 
to Eq. (1161) for the irregular input model can be con- 
structed by fixing Tf and r/, then letting t s — > oo and 
r s —5- to achieve the r; n — > limit. In this case, our 
assumption that each contact is increasingly likely to 
recover between two consecutive spikes is violated and 
Eq. (ITB1) is not valid (not pictured). Regardless, we ver- 
ify numerically that Eq. (|16[) is accurate when r m is 
decreased toward zero while keeping F in fixed (Fig.EJ). 

We now discuss the statistics of the postsynaptic re- 
sponse when the effective presynaptic spiking is much 
faster than vesicle recovery {p r r m 3> 1/t u ). In such a 
regime, incoming spikes occur much more frequently 
than recovery events and synapses becomes depleted. 
As a result, the number of vesicles released over a long 
time window is determined predominantly by the num- 
ber of recovery events in that time window and largely 
independent from the number of presynaptic spikes (Fig.[Bj3) 



Short term depression and neuronal variability 



11 



A B 




1 2 3 4 5 0.5 1 1 (sec > 

Fig. 6 Vesicle release dynamics at low and high presynaptic firing rates. A) At low presynaptic rates, vesicles are 
recovered (i.e., m(t) returns to M = 3) between presynaptic spikes. Thus, the number of vesicles released by each presynaptic 
spike is approximately an independent binomial random variable with mean p r M and variance p r M(l — M). B) At high 
presynaptic rates, vesicles are released almost immediately after they are recovered. Thus, the number of vesicles released over 
a time window of length T is approximately a Poisson random variable with mean and variance TM/t u . 



[T!?ll4"3"] . The synaptic response therefore inherits the 
Poisson statistics of the recovery events so that 

lim F x (T) = l. 

For the Poisson input model, this limit can be made 
more precise in the T — > oo limit by expanding Eq. (J8| 
in terms of the parameter a — l/(p r r- m T u ) to obtain 

F x = l-2a + 4a 2 + 0(a 3 ) (17) 

which converges to 1 as r- m r u — > oo. For the irregular 
and regular input models, we verify in Fig.[5]that F x — > 
1 when Tj n is increased while keeping F m fixed. 

The time constant, t u , at which a synapse recov- 
ers from short term depression has been measured in a 
number of experimental studies and is often found to 
be several hundred milliseconds [57U59II34U23U22U28U41] . 
Therefore, for even moderate presynaptic firing rates, 
synapses are often in a highly depleted state. As dis- 
cussed above, this promotes Poisson-like variability in 
the synaptic response. This provides one possible mech- 
anism through which irregular Poisson-like firing can be 
sustained in neuronal populations [IT] . 

4 Discussion 

We used continuous time Markov chain methods to de- 
rive the response statistics of a stochastic model of short 
term synaptic depression with three different presynap- 
tic input models. We then used this analysis to un- 
derstand how the mean presynaptic firing rate and the 



variability of presynaptic spiking interact with synaptic 
dynamics to determine the mean rate of vesicle release 
and variability in the number of vesicles released. This 
analysis revealed a number of fundamental, qualitative 
dependencies of response statistics on presynaptic spik- 
ing statistics. Some of the dependencies have been pre- 
viously noted in the literature and some have not. 

The number of vesicles released over a time win- 
dow is smaller for irregular inputs than for more reg- 
ular inputs (Figs. [T] and [2]) given the same number of 
presynaptic spikes. Thus, regular presynaptic spiking 
is more efficient at driving synapses. This mechanism 
competes with a well-known property of excitable cells: 
that they are driven more effectively by irregular, posi- 
tively correlated synaptic input currents [4*51l4"Tll3"M3"B] . 
In addition, a population of presynaptic spike trains 
drives a postsynaptic neuron more efficiently when the 
population-level activity is more irregular, for example 
due to pairwise correlations 19J. Together, these results 
suggest that a postsynaptic neuron is most efficiently 
driven by presynaptic populations that exhibit small or 
negative auto-correlations, but positive pairwise cross- 
correlations. 

Our model predicts that the de-correlating effects 
of short term depression and stochastic vesicle dynam- 
ics can produce negative temporal auto-correlations in 
the synaptic response even when presynaptic spiking 
is temporally uncorrelated or positively correlated, in 
agreement with previous studies [231TT9"] . This yields a 
response Fano factor that decreases with window size, 
as observed in some recorded data [3D] . We note, though, 
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Fig. 5 Response Fano factor as a function of presynap- 
tic firing rate for three input models. Response Fano fac- 
tors calculated over large windows for A) the irregular input 
model B) Poisson input model and C) regular input model. 
Fano factors approach 1 at high presynaptic firing rates re- 
gardless of the presynaptic Fano factor (triangle on right is 
placed at F x = 1). At low presynaptic firing rates, response 
Fano factors approach the value given in Eq. (11611 (indicated 
by triangle on left). Dotted line in (B) is from closed form ap- 
proximation in Eq. ||8j and dashed line is from the expansion 
given in Eq. (ITft . 



that some parameter choices can yield positively a cor- 
related synaptic response when presynaptic inputs are 
positively correlated [2 5) and neuronal membrane dy- 
namics can introduce positive correlations to a post- 
synaptic spiking response even when synaptic currents 
are not positively correlated in time [37] . This is con- 
sistent with several studies showing positive temporal 
correlations in recorded spike trains [^IITSIIBIITT] . 

We predict that moderate or high firing rates can in- 
duce a Poisson-like synaptic response even when presy- 
naptic inputs are non-Poisson (Fig. [5] and [IB])- This is 
because even moderate firing rates can deplete synapses 
and depleted synapses inherit the Poisson-like variabil- 
ity of synaptic vesicle recovery (Fig. [6)3 and [I9ll43j ) . 
At lower firing rates, short term depression and synap- 
tic variability can increase or decrease Fano factor. For 
example, in Fig. 03, the response Fano factor is larger 
than the presynaptic Fano factor (Fi n = 1) at low firing 
rates, decreases at higher firing rates, then approaches 
F x = 1 at higher firing rates. This complex dependence 
of firing rate on Fano factor might be related to the 
stimulus dependence of Fano factors observed in sev- 
eral cortical brain regions 

Our conclusion that fast presynaptic spiking causes 
Poisson-like variability in the synaptic response relied 
on the assumption that vesicle recovery times are ex- 
ponentially distributed. The exponential distribution 
is a justifiable choice for recovery times only if recov- 
ery times obey a memoryless property: having already 
waited t units of time for a recovery event, the proba- 
bility of waiting an additional s units of time does not 
depend on t. The precise mechanics of vesicle re-uptake 
and docking determine whether this is an appropriate 
assumption. If recovery times have a different probabil- 
ity distribution, then the synaptic response will inherit 
the properties of this distribution at high presynaptic 
firing rates instead of inheriting the Poisson-like nature 
of exponentially distributed recovery times. 

Previous methods have been developed to analyze 
the synaptic response of the model used here. In [18j 
[24] . the model restricted to the M = 1 case is ana- 
lyzed for presynaptic spike trains that are renewal pro- 
cesses. This includes the Poisson and the regular input 
model discussed here, but excludes the irregular input 
model in which the spike train is a non-renewal inho- 
mogeneous Poisson process. In [43] approximations are 
obtained for the case where the presynaptic spike train 
is an inhomogeneous Poisson process, but the approx- 
imation is only valid when the rate-modulation of the 
Poisson process is small compared to the average fir- 
ing rate. Thus, these approximations are only valid for 
the irregular input model when rt — r s <C r s . Other 
studies [32 , 36 use a deterministic synapse model that 
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implicitly treats the number of available vesicles as a 
continuous rather than a discrete quantity. This deter- 
ministic model represents the trial average of the model 
considered here and can vastly underestimate the vari- 
ability of a synaptic response [43] . 

A more detailed synapse model allows for multiple 
docking sites at a single contact [521[T5] . This model can 
yield different response properties than the model used 
here in certain parameter regimes P~5] . Even though this 
more detailed model can be represented as a continuous- 
time Markov chain, the analysis of this model would be 
significantly more complex than the analysis considered 
here since it would be necessary to keep track of the 
number of readily releasable vesicles at each contact 
separately. This would result in a Markov chain with 
K x N M states where M is the number of contacts, N 
is the number of docking sites per contact and K is the 
number of states used for the presynaptic input model 
(K = 1 for the Poisson input model, K = 2 for the 
irregular input model, and K = 9 for the regular input 
model). 

To quantify the synaptic response to a presynaptic 
spike train, we focused on the statistics of the num- 
ber of vesicles released in a time window. Postsynaptic 
neurons observe changes in synaptic conductance in re- 
sponse to presynaptic spikes. The synaptic conductance 
are often modeled in such a way that they can be easily 
derived from our process x(t) through a convolution: 
g(t) = J Q x(t — s)a(s)ds where g(t) is the synaptic con- 
ductance elicited by a presynaptic spike train and a(s) 
is a kernel representing the characteristic postsynaptic 
conductance elicited by the release of a single neuro- 
transmitter vesicle. Since this mapping is linear, the 
statistics of g(t) can easily be derived in terms of the 
statistics oi x(t) [53114*5] . 
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